eststo clear


if("`c(os)'"=="Windows"){
	local tableSaveDir = "Tables\DiD\"
	local plotSaveDir = "Plots\"
	local pythonSaveDir = "PythonScripts\"
}
else{
	local tableSaveDir = "Tables/DiD/"
	local plotSaveDir = "Plots/"
	local pythonSaveDir = "PythonScripts/"
}

//Establish Control Variable Sets - VMT
local Spec1ControlsVMT = "c.normYear#i.stateGroup c.quadYear#i.stateGroup asinhMetroPop asinhNonMetroPop asinhMetroInc asinhNonMetroInc logrealmeangasprice logemployment loglicenseddrivers logrealstategdp logroadmileage"
local Spec2ControlsVMT = "`Spec1ControlsVMT'"
local Spec3ControlsVMT = "`Spec2ControlsVMT' i.vmtDummyExtSpec3"
local Spec4ControlsVMT = "`Spec2ControlsVMT' i.vmtDummyExtSpec4"


//Establish Control Variable Sets - Gas
local Spec1ControlsGas = "c.normYear#i.stateGroup c.quadYear#i.stateGroup asinhMetroPop asinhNonMetroPop asinhMetroInc asinhNonMetroInc logrealmeangasprice logemployment loglicenseddrivers logrealstategdp logroadmileage"
local Spec2ControlsGas = "`Spec1ControlsGas'"
local Spec3ControlsGas = "`Spec2ControlsGas' i.gasUseDummyExtSpec3"
local Spec4ControlsGas = "`Spec2ControlsGas' i.gasUseDummyExtSpec4"


//Absorb Vars are always time and state fixed effects
local absorbVars = "ib51.stateGroup ib1995.year"

//Cluster variable declaration
local clusterVar = "stateGroup year"

//Treatment variable 
local treatmentVar = "1.nosafetyind"

// //Run Spec2 regressions - prep for suest command
//
// qui regress logvmt i.nosafetyind `Spec2ControlsVMT'
// scalar travelK = e(df_m)+1
// estimates store travelEq
//
//
// qui regress loghighwaygasuse i.nosafetyind `Spec2ControlsGas'
// scalar gasK = e(df_m)+1
// estimates store gasEq
//
// //Simultaneously estimate reduced form esquations
// qui suest travelEq gasEq, vce(cluster `clusterVar')
// scalar nClust = e(N_clust)
// scalar N = e(N)
// matlist e(b)["y1", "travelEq_mean:`treatmentVar'"]
// scalar delta = e(b)["y1", "travelEq_mean:`treatmentVar'"]+0
// scalar theta = e(b)["y1", "gasEq_mean:`treatmentVar'"]
// scalar deltaVarUncorrected = e(V)["travelEq_mean:`treatmentVar'", "travelEq_mean:`treatmentVar'"]
// scalar thetaVarUncorrected = e(V)["gasEq_mean:`treatmentVar'", "gasEq_mean:`treatmentVar'"]
// scalar covarUncorrected = e(V)["gasEq_mean:`treatmentVar'", "travelEq_mean:`treatmentVar'"]
//
// //Collect all necessary statistics
// scalar travelCorrection = (N-1)/(N-travelK+(nClust-1))
// scalar gasCorrection = (N-1)/(N-gasK+(nClust-1))
// scalar deltaVar = deltaVarUncorrected*travelCorrection
// scalar thetaVar = thetaVarUncorrected*gasCorrection
// scalar covarHigh = max(travelCorrection*covarUncorrected,gasCorrection*covarUncorrected)
// scalar covarLow = min(travelCorrection*covarUncorrected,gasCorrection*covarUncorrected)
// scalar covar = covarLow

//Run travel model to get appropriate N and K values.
qui reghdfe logvmt nosafetyind `Spec2ControlsVMT', absorb(`absorbVars') cluster(`clusterVar')
local k_t = e(df_a)+e(df_m)
local N_t = e(N)
local r2_t = e(r2)
local N_clust_t = e(N_clust)
local correction_t = (`N_clust_t')/(`N_clust_t'-1)*(`N_t'-1)/(`N_t'-`k_t')

//Run gas model to get appropriate N and K values.
qui reghdfe loghighwaygasuse nosafetyind `Spec2ControlsGas', absorb(`absorbVars') cluster(`clusterVar')
local k_g = e(df_a)+e(df_m)
local N_g = e(N)
local r2_g = e(r2)
local N_clust_g = e(N_clust)
local correction_g = (`N_clust_g')/(`N_clust_g'-1)*(`N_g'-1)/(`N_g'-`k_g')

//Partial urban model variables
qui reghdfe logvmt `Spec2ControlsVMT', absorb(`absorbVars') cluster(`clusterVar') res(logvmt_res)
qui reghdfe nosafetyind `Spec2ControlsVMT', absorb(`absorbVars') cluster(`clusterVar') res(nosafetyind_res_t)

//Partial urban model variables
qui reghdfe loghighwaygasuse `Spec2ControlsGas', absorb(`absorbVars') cluster(`clusterVar') res(loghighwaygasuse_res)
qui reghdfe nosafetyind `Spec2ControlsGas', absorb(`absorbVars') cluster(`clusterVar') res(nosafetyind_res_g)

//Estimate partialled urban model
regress logvmt_res nosafetyind_res_t, nocon
estimates store travel

//Estimate partialled rural model 
regress loghighwaygasuse_res nosafetyind_res_g, nocon
estimates store gas_use

suest travel gas_use, vce(cluster stateGroup)
mat var1 = e(V)
suest travel gas_use, vce(cluster year)
mat var2 = e(V)
suest travel gas_use, vce(robust)
mat mat_beta = e(b)
mat varR = e(V)

mat mat_var = var1+var2-varR

matlist mat_beta
matlist mat_var

mat diff_r = [1, 0, -1, 0]
matlist diff_r

mat diff_value = diff_r*(mat_beta')
matlist diff_value

mat diff_var = diff_r*mat_var*(diff_r')
matlist diff_var

local chi_stat = (diff_value[1,1]^2)/sqrt(diff_var[1,1])

di "`chi_stat'"

local p_value = 1-chi2(1, `chi_stat')

di "`p_value'"

// //Calculate wald-test stastics
// matrix r = [1, 0, -1, 0]
//
// matrix betaMat = [delta \ theta]
// matrix covMat = [deltaVar, covar \ covar, thetaVar]
//
// matrix outer = r * betaMat
// matrix inner = r * covMat * r'
//
// matrix chi2Val = outer*inv(inner)*outer'
//
// scalar pVal = 1-chi2(1, chi2Val[1, 1])
// //Display p-value
// di pVal
